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Abstract 

In this paper, we develop highly accurate Nystrom methods for the volume integral 
equation (VIE) of the Maxwell equation for 3-D scatters. The method is based on 
a formulation of the VIE equation where the Cauchy principal value of the dyadic 
Green’s function can be computed accurately for a finite size exclusion volume with 
some explicit corrective integrals of removable singularities. Then, an effective inter¬ 
polated quadrature formula for tensor product Gauss quadrature nodes in a cube is 
proposed to handle the hyper-singularity of integrals of the dyadic Green’s function. 
The proposed high order Nystrom VIE method is shown to have high accuracy and 
demonstrates p-convergence for computing the electromagnetic scattering of cubes 
in M. 

Key words: Electromagnetic (EM) scattering, volume integral equation, Gauchy 
principal value, dyadic Green’s function, Nystrom methods. 


1 Introduction 


Electromagnetic (EM) wave scattering in the presence of (random) mi¬ 
crostructure on material interfaces has a wide range of applications. For ex¬ 
ample, the interaction of light between the surface plasmon, a collective fluctu¬ 
ation of electron density under the background of the positive nucleus charges, 
on the metallic surfaces produces a so-called surface plasmon polariton (SPP) 
umi, which has important applications in solar cells [2], meta-materials, and 
super-resolution imaging devices W- Additionally, surface enhanced Raman 
scattering (SERS) [12] is also closely related to the excitation of surface plas- 
mons on rough or nano-pattern surfaces by incident light and is an extremely 
useful tool in hnger-printing the chemical components of a molecule, single 
molecule detector, detection of DNA studies, and bio-sensor, etc [6]. In all 
these applications, it is critical to have highly accurate and efficient numerical 
methods for computer simulations of the EM scattering in microstructures. 

In this paper, we will present a high order Nystrom volume integral equa¬ 
tion (VIE) method for the time harmonic Maxwell equations using dyadic 
Green’s functions GE(r', r). In most of the applications, the scatter is embed¬ 
ded into either a homogeneous or layered media and a Lippmann-Schwinger 
type of VIE, which is a well-conditioned second kind of Fredhohm integral 
equation, can be derived for the regions occupied by the scatter while the 
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dyadic Green’s function GE(I■^r) will ensure that the scattering held, ex¬ 
pressed in terms of equivalent current sources inside the scatter, satishes in¬ 
terface conditions along the layered interfaces as well as the Sommerfeld ra¬ 
diation conditions at oo. In the VIE formulation, the electric held inside the 
scatter will use the Cauchy Principal Value (CPV or simply p.v.) associated 
with the dyadic Green’s function as indicated in ([19]). Therefore, one of the 
most difficult issue is how to compute accurately and efficiently the CPV for 
the dyadic Green’s function which has an O(^) singularity, i.e. 

p.v. / dr' ia;Ae(r')E(r') ■ GE(r',r) = lim / dr' ia;Ae(r')E(r') • GE(r',r), 

Jo. Jn\Vs 

( 1 ) 

where Ae(r') describes the scatter dielectric constant deviating from the sur¬ 
rounding background material. As the CPV is dehned through a limiting pro¬ 
cess of diminishing size 5 of the exclusion volume, in practical computation, 
a hnite 5 has to be taken, namely, a small, hnite and hxed <5 > 0 for the 
exclusion Vs is selected, then the following approximation is taken 

p.v. / dr'ia;Ae(r')E(r') • GE(r', r) / dr'ia;Ae(r')E(r') ■ GE(r', r). (2) 
Jn Jn\Vs 

Therefore, we are faced with two issues in the implementation of the VIE: 

Firstly, what size of 6 should be taken? What is the effect of error by 
using a hnite 6 in the calculation of the CPV? For any hnite S, there will be 
a truncation error which we will call correction terms 

p.v. f dr' in;Ae(r')E(r') • GE(r',r) = f dr' in;Ae(r')E(r') ■ GE(r', r) 

Jn Jn\Vs 

+ correction terms. (3) 

The existence of these correction terms and their magnitude will limit the 
accuacy of the VIE solution if they are not explicitly included in the numerical 
solution process. The correction terms were derived by Fikioris [Tj using the 
mixed potential formulation of the electric held. In this paper, we will re-derive 
the VIE similar to those in [1], however, in a more succinct manner and the 
hnal formula is more suitable for numerical implementations. Previous work 
on how to handle the singular integrals for VIE method include singularity 
substraction [7], locally corrected Nystrom scheme [9], direct integration of 
the singularity [H], etc. 

Secondly, for the selected 6, how to hnd an accurate integration quadrature 
formula to compute the resulting integration over the domain with a 

highly singular integrand? 

It is the objective of this paper to address these two issues and hnd eas¬ 
ily implementable solutions. To address the hrst issue, we will derive the VIE 
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equation using vector and scale potentials such that the Cauchy principal value 
can be computed as in ([3]) with explicit expression for the correction terms. On 
the other hand, to address the second issue, we will construct special quadra¬ 
ture weights for the tensor product Gauss quadrature nodes in a reference 
element O (taken to be a cube in this paper) by using an interpolation ap¬ 
proach. In this approach, first a brute force computation of the integral, using 
Gauss quadrature in polar coordinate centered at the singularity, will be done 
to a given accuracy, albeit involving large number of values of the integrand. 
Then, realizing the integrand in the VIE matrix entries, except for the singu¬ 
lar denominators involving 1 < fc < 3, are smooth functions, which can 
be in fact accurately interpolated using values on the tensor product Gauss 
nodes inside the cubic domain G. Then, the brute-force integration formula 
can be converted into a new integration weights for the tensor product Gauss 
nodes. The new integration formula can be tabulated for integrating general 
functions. Finally, the Nystrom collocation method is used to discretize the 
VIE with a high order accuracy. 

The rest of the paper is organized as follows: Section 2 presents the for¬ 
mulation of a VIE where the GPV can be computed with a finite exclusion 
volume accurately. Then, numerical algorithms, including the Nystrom col¬ 
location method, efficient quadrature formula and numerical implementation 
are given in Section |3i Section 0] includes various numerical simulation tests 
such as the accuracy of Gauchy principal value computation, h-independence 
of matrix entries and p-re£nement convergence for the VIE. The paper ends 
with a conclusion in Section [5l 


2 Volume Integral equations for the Maxwell equations 


2.1 VIE and Cauchy principal values 


In this section, we will follow [3] to show briefly how the VIE for the 
Maxwell equations can be derived using a vector form second Green’s identity 
for the following vector wave equations, 

£E(r) — Ci;^e(r)E(r) = —icnJe(r), r G M^\(S U 5f2), (4) 
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where oj is the frequency, /r is the permeability, and S consists of any possible 
interfaces of the background medium in case of a layered material, 

£ = V X — Vx, 
h 

and Je(i‘) is the far-held source (assumed to be away from the layered struc¬ 
ture), which produces the incident waves impinging on the layered structure 
from the top, i.e.. 


Ei“(r) = -icn/i(r) [ Gi,(r, r') ■ Je(r')dr', (5) 

and G£;(r, r') is the dyadic Green’s function for the layered media. A scatter G 
is characterized by a different dielectric constant from the layered background 
dielectrics, i.e., 

e(r) = eL(r)-F Ae(r), (6) 

where Ae(r) = 0, r ^ G. Then, (|1]) can be rewritten as 

/:E(r) - a;^eL(r)E(r) = -icnJ(r), (7) 

where 

J(r) = Je(r) + Jeq(r), (8) 

and the equivalent current source Jeq(r) is dehned to rehect the existence of 
the scatter G: 

Jeq(r) = ia;Ae(r)E(r). (9) 


Let us consider any interior point inside the scatter, i.e., r' G G and a 
small volume Vs = Vs{y') C G centered at r'. The dyadic Green’s function 
Ge (r, r') is dehned by 


£Gij(r, r')- a;^eL(r)Gi7;(r, r') 


—-I5(r-r'), r G Ml 

/^(r) 


( 10 ) 


In the case of a homogeneous medium, we have 


GE(r.r') = GE(r',r) = (l + ivv)9(r,r'). (11) 

where = eL{r),^ is the permeability of the medium, and 


g{r, r') 


Q—ikR 

dvr R 



( 12 ) 


Next, on multiplying ([7]) by GE(r,rO and flTOl) by E(r) and forming the 
diherence, and then integrating over the domain M^\V 5 , after some manipula- 
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tion [3], we arrive at the following equation (after switching r and r'): 


iun{r) / dr' GE(r, r') ■ J(r') —/i(r) / ds' ici;GE(r, r') ■ (n x H(r')) 

Jk.^\Vs JSs 

X G£(r, r') ■ (n x E(r')) =0, r G hi, (13) 


where Ss = dVsir), n is the normal of Ss pointing out of V 5 (r). 

As (5 —?■ 0, the hrst integral will approach the Cauchy principal value 
of a singular integral, while the surface integrals will in fact depend on the 
geometric shape of the volume Vs. 

In order to estimate the surface integrals, we have the following asymp¬ 
totics for small kR 1: 


where u = (r'— r)/i?, which implies that: 

lini ^ ds' n X E(r') ■ V x Gij(r', r) = - [I - Ly,] • E(r), 

lim [ ds' n X H(r') • GE(r', r) = -—Lvg ■ V x H(r), 
< 5^-0 JSs k^ 


(14) 

(16) 

(16) 

(17) 


and the L-dyadics for Vs of various geometric shapes are given as follows |16j : 

1 . 


Lv, = 


(18) 


for a sphere. 

Substituting ffT6|) and ffT7|) into ffT3|) . after some manupulation and also 
using Ampere’s law, we have the VIE for the electric held for r G hi: 

C ■ E(r) = E’“(r) - ia;/i(r) p.v. / dr' iu;Ae(r')E(r') ■ r), (19) 

Jn 

where the coefficient matrix is given by 

C = I + Ly, ■ Ae(r). (20) 

As mentioned in the introduction, in most numerical implementations of 
([T9]), the CPV integral is computed by selecting a hnite 5 as in ([2|), therefore 
the numerical solution thus obtained does not satisfy the original VIE (IT^ 
and it will have an intrinsic error rehected in the correction terms indicated 


m 
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2.2 Reformulation of the VIE and computing CPV with a finite 6 


In this section, we will derive a VIE of the E-field where the CPV in fl20|) 
can be computed with a hnite exclusion volume together with some correction 
terms. Based on the Helmholtz decomposition, the electric held E(r) can be 
expressed as follows |3], 

E = -icuA - VV, 

according to the Lorentz gauge [13], 


V ■ A = —icue/iV, 


( 21 ) 


SO we have a vector potential representation for the electric held. 


E = —icuA 


icne/i 


-V(V ■ A) = -iw 


I + 


A. 


( 22 ) 


On the other hand, it can be shown that the potential A satishes the 
Helmholtz equation componentwise [3| 


V^A + k^A = -/iJ. 


(23) 


Thus the solution A of Eq. (12311 can be rewritten as an integral represen¬ 
tation: 


A = /X f dr'J{r')g{r,r') 

= n[ Je(r')^(r,r')dr' + /i / Jeq(r')^(r, r')dr' 
jR^\n Jn 

= jj. [ 3e{r')g{r, r')dr' + jj. [ dr'ici;Ae(r')E(r') 5 ((r, r'), (24) 

jR3\o Jn 

where the second equality on the right hand side of Eq. fl2Tll is due to the 
assumption that supp(Je(r)) fl O = 0. 

For the hrst integral in Eq. fIMll . it is well-dehned if r G O and when we 
plug it into Eq. fl22l) . it yields the incident wave E™'^(r) according to relation 
l|5|). For the second integral over O , we split it as follows 

p [ dr'icnAe(r')E(r')^(r,r') = r( f + / ) i^^Ae(r')E(r')^(r, r'), 

Jn \In\Vs Jvs J 
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and along with the hrst integral term in Eq. , it follows from Eq. (12^ that 


E = - iw 


ln\Vs 


ia;Ae(r')E(r') 


I+pVV 


r) 


lU 


I+^vv 


I Vs 


= E'^fr) - icn 


ln\Vs 


lU 




I Vs 


dr'ici;Ae(r')E(r') 5 f(r, r') 
ici;Ae(r')GE(r, r') ■ E(r') 
dr'ici;Ae(r')E(r') 5 ((r, r'). 


Next, we isolate the singular part from g{r, r') as follows 

5 '(r,r') = 5fo(r,r') +^(r,r'), 

where 


^o(r,r') = 


1 


and use the fact that m 
VV [ dr ^ 


dvr I r — r' I ’ 


= -L 


47r|r - r'| JaVs 47r|r - r' 


Vs, 


we can compute the following integral in the following manner 


(25) 


(26) 

(27) 


(28) 


VV / dr'Ae(r')E(r') 5 fo(i', r') 

Jvs 

= VV / dr'——--Ae(r)E(r) + / dr'VV 5 ^o(i', r') [Ae(r')E(r') — Ae(r)E(r)] 

Jvs 47r|r — r I Jvs 

= -Ly, A6(r)E(r) + [ dr'VVr?o(r, r') [Ae(r')E(r') - A6(r)E(r)], (29) 

JVs 


where the second integral has a removable singularity O through the 

use of the spherical coordinate transform centered at r, provided the function 
Ae(r)E(r) is a differentiable inside Vs, which we assume to be. 

Using Eq. (129|) and noting that g = g — Qq is a. regular function, Eq. (1251) 
becomes 

C . E = E'“(r) - f za;Ae(r')GE(r, r') ■ E(r') 

Jn\Vs 

+ a;^p f dr'Ae(r')E(r') 5 ((r, r') 

JVs 

2 

+ ^/i f dr'Ae(r')VV^(r, r') • E(r') 

/v JVs 
2 

+ Jy dr'VV^o(r, r') [Ae(r')E(r') - Ae(r)E(r)], (30) 

8 













with the same coefficient C being defined in Eq. fl20|) . 


The VIE in fl30p is similar to those obtained by Fikioris |1], however, our 
derivation is based on a splitting of Green’s function in fl26|) and the identity 
for the Lvg in fl28|) . A comparison study between CPV formulation flT^ and 
finite exclusion volume formulation (jdni) can be found in [15]. Now expression 
(jHnil holds for any finite 5 > 0 as long as lA C G, and all integration terms 
involved on the right hand side are well-defined provided that Ae(r)E(r) is 
Holder continuous. We can see that the last three integration terms can be 
understood as the correction terms for computing the Cauchy principal value 
with a finite-sized exclusion volume Vs- It should be noted that these integrals 
are all weakly singular integrals whose singularities can be removed by the 
simple spherical coordinate transform. In particular, we can easily estimate 
their magnitude in terms of 5. Namely, 


I [ dr'A6(r')E(r')<?(r,r')| < C,\\AeE\\J^ (31) 

JVs 

\f dr'Ae(r')VV^(r,r')-E(r')| <G2||AeE||^<5^ (32) 

JVs 

f dr'VV9„(r,r') [A£(r')E(r') - A£(r)E(r)| | < CaHAfEH, „<5, (33) 

JVs 


where Gi, G2, C 3 are constants, and || • ||oo and || ■ ||i^oo represent the L°° norms 
of a function and its first derivative, respectively. 

Remark 1 Eqs. explicitly indicate the accuracy of approximating 

the Cauchy principal value (Qj) hy the integral ^ with a finite 5 > 0, i.e., the 
truncation error is of the order 0{6). Hence the numerical solution of the VIE 
will have this 0{6) truncation error in general regardless of the integration 
quadratures used if terms in Eqs. are not included. 

However, Eor the case that Vs is a ball of radius 6 centered at r, one can 
obtain a better estimate in Eq. due to the anti-symmetry of the singular 
term VV 5 'o(r, r') in the spherical coordinates, i.e. 

I / dr'VV(?o(r,r') [Ae(r')E(r') - A 6 (r)E(r)] | < G 4 I|AeE|(34) 

JVs 
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3 Numerical Methods 


3.1 Nystrom collocation method 


Now we will use the Nystrom collocation method to solve Eq. fl30|) . First, 
the computational domain (scatter) hi is divided into N cubic elements hlj 
with length aj,z = 1,2, On each element Oj, we assign M Gauss nodes 

on which M scalar Lagrange basis functions (pij^j = 1,2,3, ....M are dehned 
and vanish outside Oi. Then, we can write the solution as 


N M 

E(r) Gi0q(r), r eQi, 

i=i j=i 


(36) 


with Cij are the unknown vectorial coefficients. Inserting Eq. fl55]) into Eq. 
), we obtain the following equations for Cij-. 


N M 


c.c„=Ej“+i^VEE 


n=l m=l 




dr'Ae(r') Ge (rq-, r')(j)nm (r') 


M 


+ uj^ia 


m=l 


'Vs-' 


O M 
(jJ 11 X—^ 

+ T:rE 


m=l L' 


dr'Ae(r')g'(ry,r')0im(r') 

[ dr'Ae{r')Wg{rij,r')(l)im{r') 

JVx.. 


2 M 
ui y, 


+ I dr'V^5ro(rii,r') [Ae(r')0im(r') - Aeij(j)im{rij)] 


m=l 


(36) 


When n ^ i, we have Vs^. ^ O, then in Eq. fl36|) the integral 


InAVs, 


dr'A€(r')GE(rij, r')')>nm(r') 


(37) 


is regular on the whole cube and hence it can be evaluated by the regular 
Gauss quadratures, i.e.. 


r _ /a ^ 

/ dr'Ae(r')GE(ry,r')0„m(r') = I ^ j ^ Ae„mGE(ry, r„^)a;^, (38) 

\ / / 


with 00 ^ being the standard Gauss weights in 3-D, which are obtained from 
the tensor product of the Gauss weights in the 1-D domain [—1,1]. 


When n = i, although the singularity Cj, is excluded from the domain Dj, 
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the calculation of the integral 


[ dr'Ae{r')GE{rij,r')(j)im{r') (39) 

is still challenging. Next, we present an efficient quadrature formula to evaluate 
this integral. 


3.2 Interpolated weights on Gauss nodes for integrals on 


For the sake of generality, we consider the following integral: 

/(r;rXr;r') 


Is = 


ln\Vs 




-dr', r e Vs, 


(40) 


where k = 1,2,3 corresponds to weak, strong, and hyper-singularity of the 
integral, respectively. The function /(r; r') is assumed to be a general smooth 
and well-defined function, while h.(r; r') is a fixed arbitrary function which 
results from the directional derivative in the definition of the dyadic Green’s 
functions. 


As the function /(r;r') is smooth over the whole domain G, then it can 
be well approximated by the following simple interpolation: 

/(r; r') ^ ^ /(r; rj)0j(r'), G G, (41) 

i=i 

where are J-nodes in G formed by the tensor product of Gauss node 

in [—1,1]. Combining Eqs. fl40|) and fl4T|) yields 



/(r;r')h(r;r') 


dV 


j 

E 

i=i 


/(r;G 




(42) 


We call 00 j the effective interpolated weights, or just interpolated weights, which 
are defined through the integral 


S' = 


i 


S(r')h(r;r') 


n\Vs 


r — r 


/1 k 


dr'. 


(43) 


Note that the interpolated weights ooj depend on the location of the singularity 
r and relies on accurate calculations of Eq. (l43i) . which will be accomplished in 
two steps: first the domain G is subdivided into nine cubes with one containing 
the singularity r in its center. For the cube including the singularity, a straight¬ 
forward, brute-force approach involving a large number N {N ^ J) of Gauss 
points is adopted in local spherical polar coordinates, will be used to obtain 


11 





satisfactory accuracy. While for the other cubes, regular tensor product Gauss 
quadrature is applied. Details of the computations and resulting weight tables 
could be found in HZ]. 

Fortunately, the computation of weights only needs to be performed 
once and tabulated for the reference domain, then they can be used for varies 
cubic elements. Due to the smoothness of function /(r,r'), the number J 
is relatively small, especially if the element size is small as in practice, so 
computation of Eq. fl5U]) is now efficient once are obtained. 

Remark 2 As the Cauchy principal value in the VIE, thus the VIE itself, 
depends on the specific shape of the exclusion volume, the interpolated weights 
can only he used for the shape for which it was calculated, namely, a cubic 
shape here. So, for a general element obtained by a affine mapping as in a 
finite element triangulation, we will need to isolate the singularity by a cube 
with the singularity at its center, then the pre-calculated interpolated weights 
defined in can be used while integral over the rest of the region outside 
the cube within the element of a more general shape can be done with regular 
Gauss guadratures. 


3.3 Computation of VIE matrix entries 


In this section, we will show how to compute the matrix entries accurately 
for the VIE in the following steps. 

• Step I: calculate interpolated weights on the reference cubic domain. 


For Eq. fl43|) . we take D = [—1,1]^ and Vs = B{rj,6), where 5 > 0 is a 
prescribed small number and rj,j = 1 , 2 ,...M coincident with the coordi¬ 
nates of M Gauss points in the reference domain, and we take J = M in 
Eq. (jH]). 

Further, the dyadic Green’s function has the form 


Ge — gl + 




-,—ikR 


A-nR 


le 


—ikR 


AnR^k 


(I — 3u 0 u) — 


u 0 uj 

^—ikR 

ATiR^k^ 


l-3u( 


u 


(44) 


Since the value of function u 0 u is multiple-defined at i? = 0, in fl4U]) we 
need to take 

h(r; r') = u 0 u, (45) 

and hence Eq. fl43|) yields a set of 9 interpolated weight matrices. However, 
due to the symmetry of the matrix u 0 u, only 6 components need to be 
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calculated. For the identity matrix I term in Eq. (jHj), we will also need a 
set of scalar interpolation weights by assuming h(r, r') = 1 in fl43|) . 

Additionally, for the scalar and matrix weights, we need to calculate for 
k = 1,2, and 3 for weak, strong, and hyper singular integrals, respectively. 

To summarize, for each collocation point (also the singularity location) 
rj,j = 1,2,...,M in an element, scalar weights and = 

1 ,2,...,M are calculated for weak, strong, hyper-singular integrals, respec¬ 
tively. And the corresponding weights are denoted as A^ m, ^j,m, and Aj^m- 
These weights only need to be calculated once and then stored for future 
use. 


• Step II: We discretize the computational domain into a series of cubic el¬ 
ements with length ai,i = 1,2, ...,iV and assign M Gauss points in each 
elements. For the j-th Gauss point rjj in the i-th element, we construct the 
equation: 

N M M / 1 \ 

^ h ^ ^ ^ ) -^nm ■ C-nm ^ ) ^im ' ^im ^ q ^*^*1 j ^3x3 ' ^ij ^ij ’ 
n=l m=l m=l \ O / 

The matrix B originates from the correction terms of the Gauchy principal 
value 


B 


im 


ui'^H [ dr'Ae{r')g{rij,r')(j)im{r') 

^ [ dT'Ae{r')V^g{rij,T')(pim{r') 

2 

^ [ dr'V^5fo(rp,r') [Ae(r')0im(r') - Aeijcpimirij ))], (47) 

K JB{rij,ai5) 


and it can be calculated by standard Gauss quadrature through spherical 
coordinates since the Jacobian will eliminate completely the singularity of 
the integrands. 

When n = i, we calculate the integral Eq. fl39l) as 



where Rm = |rp — fiml and recall the definition of the interpolated weights 
Eq. fl4^ . we have 


= 

A* = 


^j,m, 


ii) 

(^]a- a* =(^Ya- a* 


ay 


^j,mi ^j,m 
2 - 



(49) 
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when n ^ i, we have 


A 


nm 



(50) 


• Step III: Eq. fH6l) for all the N x M Gauss points can be assembled as the 
following linear algebraic equation system 



^ XX 

^xy 

V,, 




ET 

V c = 

^yx 

^yy 

V,. 


Cy 

= 

Einc 


V 

^ zx 

V 

v zy 

V,, 


Cz 


Einc 


Based on the properties of the Green’s function, the 3NM x 3NM matrix 
V is partitioned into nine blocks, each of which is a. NM x NM sub-matrix. 
The solution of the VIE contains three iVM x 1 vectors, which represents 
the held in x, y, and z directions. This system is solved by iteration methods 
such as the GMRES method. 


4 Numerical Results 


In this section we test the accuracy of the interpolated weights on Gauss 
nodes, the (5-independence of the solution of the VIE, and convergence of the 
p-rehnement of the Nystrom collocation method. 


4-1 Accuracy of the interpolated weights on Gauss nodes 


In Eq. fl46|) . the calculation of matrix B from the correction terms are 
straightforward, so we will focus on validating the interpolated weights in com¬ 
puting matrix A. For convenience, we consider the integral of a real-valued, 
tensor function 


cos R 
R 


u (g) u) 


cosi? 

R2 


3u (g) u) -f 


cos R 
R3 


3u(g) u). 


(52) 


which is similar to the Green’s function Eq. fl44|) on the domain G\V 5 . Without 
loss of generality, we take G = [—1,1]^ and rj as the 27 points constructed 
from the tensor product of the Gauss points of order 3 in [—1,1]. Thus, we 
have j = m = 1, 2,..., 27 as in Eq.(jl9]) and the integral yields 


27 


m=l 


^ ^ ^ COsdr^TT, “1“ ^j,m) I ; (^3) 
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Table 1 

Convergence of the integral as the singularity is at center, gn = §22 = <733 and 
512 = 513 = 523 = 0- _ 



6 = 0.1 

(5 = 0.05 

(5 = 0.025 

6 = 0.0125 

5ii 

3.985701 

4.017024 

4.024872 

4.026835 

error 

-4.1784E-2 

-1.0461E-2 

-2.613E-3 

-6.5E-4 

order 

- 

2 

2 

2 


Table 2 

Convergence of the integral as singularity is at a corner. Reference solution gn = 
522 = 533 = 0.982526 and gi 2 = 513 = 523 = —0.998097. _ 



5 = 0.1 

6 = 0.05 

5 = 0.025 

6 = 0.0125 

5ii 

0.940714 

0.972063 

0.979913 

0.981876 

error 

-2.80425E-1 

-7.3424E-2 

-1.8102E-2 

-3.983E-3 

order 

- 

1.93 

2 

2.1 

512 

-0.998084 

-0.998094 

-0.998097 

-0.998097 

error 

1.3E-5 

3.0E-6 

0 

0 

order 

- 

2.1 

- 

- 


which is a 3 X 3 matrix depending on r^. 

For each Gj, we use the direct brute force method introduced in ca to 
obtain the reference solution with the very small 6 = 10“^. Then we calculate 
the integral using the interpolated weights as in Eq. flS^ with different values 
of 6 . According to the previous analysis, the differences with the reference 
solution should decay in the order 0 {6‘^) as 6 decreases, which will be checked 
in the following. 

We classify the 27 sets of weights into four categories, based on the position 
of the singularity r^, as the ones near the corner, edge, face, and center of the 
cube. The matrix Gj is symmetric, so we only check the three diagonal entries 
idn, 922, and 533) and the three upper diagonal entries (512, 513, and 523)- 

Table [T] presents the numerical results when the singularity rj is located 
in the center of the cube, in which case 511 = 522 = 533 and the off-diagonal 
entries are all zeros. The values of 6 are taken as 0.1, 0.05,0.025, and 0.0125 
while the reference solution is 511 = 4.027477. 

In a similar fashion. Tables [2ll4] show the accuracies when the singularity 
is located near the corner, edge, and face of the cube, respectively. Comparing 
to the reference solutions, the expected second-order decay with respect to the 
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Table 3 

Convergence of the integral as singularity is at an edge. Reference solutions gn = 
-1.515302, ff 22 = g 33 = 3.39234, ff 23 = -1.579086 and gi 2 = gi 3 = 0- 



6 = 0.1 

(i = 0.05 

(i = 0.025 

(5 = 0.0125 

gu 

-1.559532 

-1.526351 

-1.518059 

-1.515987 

error 

-4.423E-2 

1.1049E-2 

2.757E-3 

-6.85E-4 

order 

- 

2 

2 

2 

922 

3.35175 

3.38217 

3.389798 

3.391707 

error 

-4.059E-2 

-1.1017E-2 

2.542E-3 

-6.33E-4 

order 

- 

1.9 

2.1 

2 

923 

-1.579072 

-1.579082 

-1.579085 

-1.579085 

error 

1.4E-5 

4.0E-6 

l.OE-6 

l.OE-6 

order 

- 

1.8 

2 

0 


Table 4 

Convergence of the integral as singularity is at a face. Reference solutions gn = 
0.877428, g 22 = 0, g 33 = 6.494784, and gi 2 = giz = g 23 = 0- _ 



S = 0.1 

6 = 0.05 

6 = 0.025 

6 = 0.0125 

511 

0.83442 

0.866672 

0.874742 

0.87676 

error 

-4.3008E-2 

-1.10756E-2 

-2.686E-3 

-6.68E-4 

order 

- 

1.9 

2 

2 

933 

6.455419 

6.484909 

6.492315 

6.49417 

error 

-3.9365E-2 

-9.875E-3 

-2.469E-3 

-6.14E-4 

order 

- 

2 

2 

2 


6 is confirmed. 


4-2 6-independence of the VIE solution 

Equation (130|) provides a formulation from which the solution of the VIE 
will be independent of the choice of a specihc 6 . In the following tests, we take 
p = 1, Ae = 4, a; = 1 and solve the VIE. The computation domain is taken 
as [—7r/2,7r/2]^, while the incident wave is 

■pine _ ik(—y-\-0.5z) pine _ pine _ n 
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Column number in a row of matrix entries Column number in a row of matrix entries 

(a) (b) 

Fig. 1. Differences of matrix entries with 5=0.1 and 5 = 0.001. (a): without correc¬ 
tion terms; (b) with correction terms. 

We first check the ^-independence of the matrix entries in Eq. (15T]) . Figure 
[T] displays the differences of one row of entries in the matrix V between the 
choices of 5 = 0.1 and 6 = 0.001, in which the solid lines are for the entries from 
a diagonal block ( 14 a;) and dashed lines are for the entries from an off-diagonal 
block {Vxy). The blue curves are for real while red curves are for imaginary 
parts. From Fig. [T](a) we can see that the differences between entries in the 
corresponding positions can be as large as 3.0 x 10“^ when the correction 
terms are not included. In contrast, the corresponding differences are reduced 
to below 2 X 10“^^ when the corrections are, hence the matrix entries are 
h-independent. 

Next we check the h-dependence of the overall solution of the VIE. The 
solution of VIE with a very small 6 = 0.001 is taken as the reference solution. 
Then various choices of 6 are taken and the differences of the corresponding 
solutions with the reference solution are calculated. The differences are mea¬ 
sured in the L°° norm for the three components Ex, Ey, E^ and they are listed 
in Tables [SMl 

Table 5 

Comparison of solutions of the VIE without the correction terms 



5 = 0.1 

5 = 0.05 

6 = 0.025 

5 = 0.0125 

WEx-EfU^ 

WEy-EfU^ 

II rp ZTTef II 

3.360E-3 

1.476E-3 

2.533E-3 

8.264E-4 

3.696E-4 

6.387E-4 

2.044E-4 

9.258E-5 

1.597E-4 

5.039E-5 

2.358E-5 

3.969E-5 


From Table [5] it can be seen that without the correction terms, the solution 
of VIE has an obvious dependence on the choice of 5 and the differences 
follow the order of 0 {6^), while the solution is indeed h-independent when the 
correction terms are included, as shown in Table El 
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Table 6 

Comparison of solutions of the VIE with the correction terms 



6 = 0.1 

(5 = 0.05 

6 = 0.025 

6 = 0.0125 

W^X - 

TTTOf 11 

- \\L^ 

8.0E-12 

l.OE-12 

0 

0 

W^y ~ 

Tprei II 

-Ey ||L°° 

2.0E-12 

l.OE-12 

0 

0 

\\Ez- 


l.OE-12 

0 

0 

0 


4-3 p-convergence of the VIE 


For numerical convergence, one can discretize the computational domain 
into finer elements (increasing the number N, h-refinement) or use higher order 
polynomial basis (increasing the number M, p-refinement) for the integral in 
the VIE. In the current work we focus on the latter since the emphasis is on 
the accurate calculation of the Cauchy principal value of the integral with 
singularities in the domain. 

In the following tests, the Lagrange interpolation of the solution of the VIE 
with M = 7 (namely, 6th order polynomial basis functions) is taken as the 
reference solution and differences between the solutions E^, M = 3,4, 5, 6 
are taken. We consider the linear relation between the logj^Q of the energy error 

Error = ||E“ - 

and the order p = M — 1 of polynomial basis function, i.e. 

log (Error) = /cp + 6, 

where k and b are parameters. In Fig. |2] there plots the logj^Q error against the 
order p, as well as the fitted lines, for the solution components E^, Ey, and E^ 
respectively. It can be seen that the log(Error) decays linearly with respect 
to p. Quantitatively, we obtain 

k = -0.464, b = -0.092, for E^ 

k = -0.353, b = -0.763, for Ey • 

k = —0.391, b = —0.946, for E^. 

Hence we conclude that the exponential convergence is achieved for the errors 
against the order p. 

Figure [3] is a 3D plot of the incident wave and resulting electric fields inside 
the cubic scatter. Rows 1-3 of Fig. [3] are for the computations with Gauss node 
M = 3,5, 7, respectively, and from the left to the right are the electric field 
components E^, Ey, E^ and the incident wave E™^, respectively. 
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(b) 


Fig. 2. logio errors of solutions agains order of polynomial basis funtions for (a), 
Ey (b) , and E^ (c). 

Figure m displays the electric field in the free space where nine cubic scat¬ 
ters are embedded. In order to clearly illustrate the wave interactions among 
the scatter array, only the scattering held E(r) — E‘“'^(r) are shown in contours 
in a specihc cross-section plane. In these tests, the incident wave is taken as 

■pine _ pine _ pine _ ik{—2x-\-2y) 

Each cube has a length of 0.5 and they form a 3 x 3 array align in the x-y 
plane. The coordinate of the center of the hrst cube is (0.25,0.25,0.25). We 
consider two cases: in the hrst case the cubes are relatively far way from each 
other (0.25 apart) and in the second case the cubes are relatively closer (0.1 
apart). In each case a lower permittivity (Ae = 4) and a higher permittivity 
(Ae = 16) are taken. Contours of the electric held in the plane z = —0.1 
are shown in Fig 01 The hrst and second rows are for the scatters with smaller 
and bigger distances, respectively, in which the left is for Ae = 4 and the right 
is for Ae = 16. It can be seen that the electric waves have stronger interactions 
when scatters are closer to each other. 


5 Conclusion 


In this paper, we have developed an highly accurate and efficient method to 
solve numerically the volume integral equation (VIE) of the Maxwell equation 
involving the Cauchy principal value of the singular dyadic Green’s function. 
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Fig. 3. Electric field and the incident wave with M = 3 (first row), M = 5 (second 
row), and M = 7 (third row). 

The VIE used allows us to compute the CPV for hnite size exclusion volume 
Vs without truncation errors (of the order of 0(5) in general and it is 0(5^) 
when the excluded volume is a sphere), which were accurately accounted for by 
including corrections terms. With the correction terms, the numerical solution 
of the VIE will be ^-independent. In addition, an efficient quadrature formula 
has been introduced to accurately compute the type integration on the 
domain f2\V5. Finally, a Nystrom collocation method based on the proposed 
quadrature formula were applied to the VIE resulting in an accurate and <5- 
independent solution of the electric held. 

The developed algorithms are verihed with several numerical tests. First, 
the accuracy of the interpolated quadrature weights for dyadic Green’s func¬ 
tion was conhrmed for various locations of the spherical exclusion volume Vs 
inside a cubic scatter. Secondly, the procedure of computing the Cauchy prin¬ 
cipal values for the dyadic Green’s function with a hnite size 6 is shown to be 
accurate with the help of the proposed correction terms. As a result, the so- 
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Fig. 4. Contours of the electric field from scatter arrays, (a)-(b): scatters are 0.1 
apart from each other; (c)-(d): scatters are 0.25 apart from each other, (a, c): Ae = 4; 
(b, d): Ae = 16. 


lution of the VIE is shown to be ^-independent, and hnally, we demonstrated 
the convergence of the VIE for p-re£nement with increasing order of basis 
functions p = m — 1, where the number of Gauss node along each direction 
m = 3,4, 5, 6, and 7, respectively. 
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